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ABSTRACT 

We use 500 pc resolution cosmological simulations of a Virgo-like galaxy cluster to 
study the properties of the brightest cluster galaxy (BCG) that forms at the center of 
the halo. We compared two simulations; one incorporating only supernovae feedback 
and a second that also includes prescriptions for black hole growth and the resulting 
AGN feedback from gas accretion. As previous work has shown, with supernovae 
feedback alone we are unable to reproduce any of the observed properties of massive 
cluster ellipticals. The resulting BCG is rotating quickly, has a high Sersic index, 
a strong mass excess in the center and a total central density profile falling more 
steeply than isothermal. Furthermore, it is far too efficient at converting most of the 
available baryons into stars which is strongly constrained by abundance matching. 
With a treatment of black hole dynamics and AGN feedback the BCG properties are 
in good agreement with data: they rotate slowly, have a cored surface density profile, 
a flat or rising velocity dispersion profile and a low stellar mass fraction. The AGN 
provides a new mechanism to create cores in luminous elliptical galaxies; the core 
expands due to the combined effects of heating from dynamical friction of sinking 
massive black holes and AGN feedback that ejects gaseous material from the central 
regions. 

Key words: black hole physics - cosmology: theory, large-scale structure of Universe 
- galaxies: formation, clusters: general - methods: numerical 
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1 INTRODUCTION 

In the standard ACDM cosmological scenario, galaxy for- 
mation proceeds hierarchically and galaxy clusters are the 
largest and most recently assembled bound systems in the 
universe. The entire Hubble sequence can be found in galaxy 
clusters; from blue extended spirals to red massive elliptical 
spheroids. For this reason galaxy clusters can be considered 
as ideal environments for the study of galaxy formation and 
evolution. The brightest cluster galaxy (BCG), is usually 
the largest galaxy that lives at the bottom of the cluster 
dark matter halo. Its stars can get there in a number of 
ways. Some of the stars will form within the rarest peaks 
collapsing at high redshift - those galaxies and their dark 
matter halos quickly virialise through violent rapid merging 
and end up at the centre of the cluster. Stars can also be 
aquired after the cluster formed through the dynamical fric- 
tion and merging with other BCGs or with the most massive 
satellites of the cluster. Finally, the stars may form in-situ 
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from the cold gas that is expected to collect at the centre of 
the potential. The satellite galaxies themselves can undergo 
further evolution inside the cluster due to the interaction 
with the hot intracluster medium, or gravitational interac- 
tions with the cluster potential and other massive galaxies. 
There's growing evidenc e that present da y massive galaxies 
are formed in two phase l|Oser et al.ll20ld h at redshift z > 2 
stars are formed within galaxies from infalling cold gas, at 
redshift z < 3 galaxies mainly grow through accretion of 
stellar material. Massive early type galaxies can assemble 
more than half of their present day mas s through dry mi- 
nor mergers at z < 2 jNaab et all 120091 ). These processes 
affect their star-formation histories and their morphological 
appearance. 

Explaining the formation of very massive BCGs at 
the center of galaxy clusters is one of the most challeng- 
ing open problems in galaxy formation research. Standard 
galaxy formation models in the context of ACDM cosmol- 
ogy are affected by the so-called "overcooling problem": 
massive galaxies are predicted to be too bright and too 
blue when compared to massive galaxies in the nearby 
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universe (|Borgani fc Kravtsovl [20091) . These models find a 
stellar content in massive cluster galaxies that is signif- 
icantly above the observed values i Kravtsov et al.l l2005h . 
even if extreme sup ernovae feedback recipes are included 
l|Borgani et al.l l2004). One scenario that has been proposed 
to solve this problem involves feedback processes from Super 
Massive Black Holes (SMBHs), usua lly referred to as AGN 
feedback. Theoretical considerations (TTabor fe Binnevll 19931 ; 
ICiotti fc Ostrikerl 19971 : ISilk fc Reesl 19981 ) suggest that these 
processes should provide enough energy to prevent gas from 
accumulating in the central regions of galaxy clusters thus 
quenching star formation. Semi-analytical models including 
AGN feedback coupled with N-body simulations have shown 
that most of the stellar mass in present time BCGs was as- 
sembled through dry minor mergers, following a phase of 
quiescent star forma tion influenced by feedback processes 
l|Lucia fc Blaizotl 120071 ). The strongest evidence supporting 
the existence of AGN feedback is provided by observations 
of X-ray cavities and radio cavities in galaxy clusters. These 
have been interpreted as buoyantly rising bubbles of high 
entropy material injected in the cluster core by jets of rela- 
tivistic particles. The success of numerical hydrodynamical 
simulations including AGN driven buoyantly rising bubbles 
in reproducing masses and colours of o bserved BCGs has 
been shown bv lSiiacki fc Springell l|2006l ). Thus, AGN feed- 
back is expected to play a significative role in shaping the 
properties of BCGs. 

In this paper we use two high resolution cosmological 
simulations of a Virgo- like cluster, with and without AGN 
feedba ck. The simulations were performed bvlTevssier et all 
d2010l) using t he A MR code RAMSES jTevssieij |2002|) . 
iTevssier et all (|2010l ) showed that combining high spatial 
resolution and the effects of AGN feedback it is possible to 
bring the stellar content in the cluster closer to the observed 
values, while a model without AGN feedback totally fails. 
Here, we use the same simulations to study the effect of 
AGN feedback on the massive BCGs that form at the cen- 
ter of the cluster. The main questions we want to address 
are: 

• What are the differences between the properties of the 
BCG when we include AGN feedback and when we do not? 
We consider structural as well as kinematic quantities to 
address this problem, including ellipticities, masses, velocity 
dispersions, rotational velocities and stellar surface density 
profiles. 

• What are the main ingredients that determine the 
evolution of the BCG? Massive elliptical galaxies are ex- 
pected to form through a long series of dry minor mergers 
i|Naab et al]|200ll2009l ). but processes connected to SMBHs 
and AGN feedback may also influence their evolution. 

• Recent results showed that it is possible to repro- 
duce some of the properties of massive elliptical galaxies 
in the field and in galaxy groups neglec ting AGN feedback 
l|Naab et al.ll2009l : lFeldmann et alj|2010f ). Is it possible to do 
the same for BCGs in massive galaxy clusters? If not, can 
the problem be solved by invoking AGN feedback? 

The paper is organized as follows: the first section is 
dedicated to the numerical methods and the sub-grid recipes 
adopted for our simulations (cooling, star formation and 
AGN feedback), while the second section presents our re- 
sults, comparing our models with observational data and 



with other numerical simulation studies. The final section is 
left for discussion. 



2 NUMERICAL TECHNIQUES 

In this section we describe the numerical methods and the 
initial conditions used to model the Virgo-like cluster we will 
study in this paper. We consider two numer ical cosmologi- 
cal sim ulations that have been presented bv lTevssier et all 
(2010). They were performed using the zoom-in technique 
which allows us to obtain the required effective resolution 
in selected regions of the computational domain. We con- 
sider a periodic cubic box of the universe of side 100 Mpc/h 
that has the standard ACDM cosmological parameters, with 
Sl m = 0.3, Ov = 0.7, Sl b = 0.045, ag = 0.77 and H = 70 
km/s/Mpc. We used the lEisenstein fc Hul (I1998T) tr ansfer 
function and the grafic package 1 Bertschingei I2OO1I) in its 



parallel implementation mpgrafic ( Prunet et al. 20081 ) to 



generate our initial conditions. We first ran a low resolu- 
tion dark matter only simulation, then we identified dark 
matter halos at z = 0. From the original halo catalog we 
constructed a set of candidate halos whose virial masses lie 
in the range 10 14 to 2 x 10 14 M /h. Finally, we identified 
our final halo based on its assembly history: most of its mass 
is already in place at z = 1, therefore it can be considered 
as relaxed at z = 0. Its final virial mass is M vir ~ 10 14 M , 
and A'/ 2 ooc = 1.04 x 10 14 M Q or M 50 oc = 7.80 x 10 13 M , 
where indice c refers to the critical density. We stress that 
our aim is not to exactly reproduce the properties of the 
Virgo cluster, that is a relatively unrelaxed cluster, but to 
study the effect of AGN feedback on the formation of the 
central galaxy in a cluster with a similar mass. We select a 
relatively relaxed halo: this allows us to focus on this effect 
without having to deal with a very complex merger history. 
The Virgo-like halo was then re-simulated at higher resolu- 
tion including dark matter and baryons. 



2.1 The simulations 

The two cosmological hydrodynamical simulations consid- 
ered in th is paper were performed using the AMR code 
RAMSES l|Tevssierj|2002l ). The initial grid had an effective 
size of 2048 3 and it has been used to extract a set of high 
resolution dark matter particles only in the Lagrangian vol- 
ume of the halo, while a lower resolution has been used to 
sample the rest of the periodic box. As a result we have 
22 x 10 6 particles in the cosmological box, 19 x 10 6 particles 
in the high resolution region and 8 x 10 6 particles within the 
virial radius of the halo at z = 0. This means that in the 
high resolution region we have a dark matter particle mass 
of 8.2 x 10 6 M© /h and a resolution element mass of 1.4 x 10 6 
M0/I1 for the baryonic component. 

The AMR grid used to solve the hydrodynamical equa- 
tions was initially refined to the same level as the particle 
grid (2048 3 , level I = 11), but 7 more levels of refinement 
were considered during the run (level ^ ma x = 18). We used 
a refinement criterion that allowed spatial resolution to be 
nearly constant in physical units; in this way the minimum 
cell physical size was always close to Aa; m i n = L/2* max ~ 
500 pc/h. The grid was dynamically refined using a quasi- 
Lagrangian strategy: when the dark matter or baryonic mass 
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Mass and spatial resolution 



^cdm 


'^gas 






[10 6 M /h] 


[10 6 M /h] 


[10 6 M /h] 


[kpc/h] 


8.2 


1.4 


0.3 


0.38 



Table 1. Mass resolution for dark matter particles, gas cells and 
star particles, and spatial resolution (in physical units) for our 
two simulation. 



in a cell reaches 8 times the initial mass resolution, it is split 
into 8 children cells. 



2.2 Modelling the baryonic processes 

Gas dynamics is modeled using a second-order u nsplit 
Godunov scheme i|Tevssierl 12002; iTevssier et all 120061 ; 
iFromang et al.l [20061 ) based o n the HLLC Rie mann solver 
and the MinMod slope limiter (|Toro et al.ll 1994 ) . We assume 
a perfect gas equation of state (EOS) with 7 = 5/3. Part of 
the galaxy formation process has been calculated using sub- 
grid models, since a correct modeling of the turbulent and 
multiphase interstellar medium (ISM) is beyond the reach 
of present-day cluster simulations. In both simulations: 

• The ISM has been modeled using a very simple EOS 
for the gas 



Tfloor — T* I — 



(1) 



where n t = 0.1 H/cc is the density threshold that defines 
the star forming gas, T* = 10 4 K is a typical temperature 
mimicking both thermal and turbulent motions in the ISM 
and F = 5/3 is the polytropic index controlling the stiffness 
of the EOS. Gas cannot cool below the temperature floor, 
while it can be heated above. 

• Gas cooling is followed according to the 

ISutherland fc Dopital (|l993l ) cooling function. We take 
into account H, He and metal cooling. Gas metallicity 
is advected as a passive scalar, and is self-consistently 
accounted for in the cooling function. We also considered 
th e effect of the st a ndard homogeneous UV background 
of lHaardt fc Madaul j 19961 ). but we modified the starting 
redshift, extrapolating the average intensity from z rc ion = 6 
to Zrcion = 12. This extrapolation is justified by the fact 
that early reioni zation is expected in proto-cluster regions 
jlliev et alj|2008h . 

• Star formation is implemented using the following sim- 
ple model. We create new star particles in cells with gas den- 
sity larger than n, . The mass of the star particles depends 
on resolution; in the present case we have chosen 3 x 10 5 
Mq /h. The formation rate of star particles is given by 



P* 



Pgas 
' t S 



with iff 
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(2) 



where tg is the local free-fall time of the gaseous component 
and e* = 0.01 is the star formation efficiency. 

• Super-novae feedback is implemented in the code. A 10% 
mass fraction of each star particle is ejected in supernovae 
explosions after 10 Myr. We assume that the supernova en- 
ergy is 10 51 erg. We also include metal enrichment from 



supernovae: one Mq of metals per 10 Mq average progen- 
itor mass is ejected in the ISM. This value for the metal 
yield has been chosen to ma tch, on average, the meta l en- 
richment from massive stars l|Wooslev fc WeaveJI 19951 ). Su- 
pernovae feedback_jsinodeled using the "delayed cooling" 
scheme (|Stinson et al.l [20061 ). i.e. we shut down gas cooling 
for 50 Myr in the cells surrounding a star particle going 
supernova. Since regions where feedback acts are typically 
very dense, this technique pre vents gas from fastly radiating 
away the feedback energy. See lAgertz et all l|201ll ) for more 
details. 

It has been shown that these galaxy formation recipes 
are able to successfull y reproduce the properties of spi- 
ral galaxies in the field dMaver et a l. 2008; Governato et al.l 
120091 l2010l : lAgertz et all l201ll ) as well as other observed 
galaxy properties like the Ken nicutt-Schmidt law , star 
formation rates, galactic winds jDubois fc Tevssierl 120081 : 
iDevriendt et alJl201Cl : lAgertz et al.ll201ll ). On galaxy group 
scales the same recip es are less successful in reproducing the 
observed properties ( Feldmann et a l. 2010), while on galaxy 
clust e r scales they fail (|Borgani et al.ll2004l : iKravtsov et all 
I2OO5I : iBorgani fc Kravtsovll2009l ). In simulated galaxy clus- 
ters similar sets of phenomenological models produce over- 
cooling of gas that leads to the formatio n of a higher fraction 
of sta rs than observed in real systems |Borgani fc Kravtsovl 
2009). It has been argued that star formatio n is quenched in 
halos more massive than M c ~ 6 x 10 11 Mq (|Cattaneo et al.l 
2006), and this mass threshold is thought to be related to 
the stabilisation of gas accretio n shocks and to the trans ition 
from cold to hot gas accretion dBirnboim fc Dckel 200^), al- 
though additional physical processes are required to prevent 
overcooling. The favoured theoretical scenario proposed to 
solve the overcooling problem involves the role of AGN feed- 
back in halos of mass M > M c . Gas accreti on favours AGN 
activ ity and the related feedback on the gas (|Cattaneo et al.l 
2006). AGN feedback is expected to provide enough energy 
to heat up gas in halos and partially blow it away, thus pre- 
venting further star formation. The analysis of cosmological 
simulations of galaxy g roups and clusters in cluding AGN 
feedback like the ones of lPuchwein et all (|2008l ) support this 
scenario. 

The two simulations of the Virgo-like cluster we are 
considering differ substantially. The first run has been per- 
formed using only the galaxy formation recipes described 
above, without considering the presence of SMBHs and ne- 
glecting AGN feedback, therefore after much thought of a 
clever acronymn, we call it the AGN-OFF run. In the sec- 
ond run we include SMBHs and we take into account AGN 
feedback; we call it the AGN-ON run. 



2.3 SMBH growth and AGN feedback in the 
simulations 

The seeds for SMBHs formation are thought to be either Pop 
III stars (|Madau fc Roes 2001), or a result direct collapse 
of baryonic material within low angular mo mentum halos 
(Bro mm fc Loeb 2003; Bcgcl man et al.ll2006l ). In both cases 
the seed SMBHs are expected to grow relatively quickly to 
Mbh, s = 10 5 M when they will start to interact with the 
environment and self-regulate their gas accretion rate. The 
fact that this black hole mass is at least one order of mag- 
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nitude lower than th e minimum SMBH mass observed in 
the M B h - cr relation l|Gebhardt et alj|200d : iGiiltekin et all 
2009) allows us to consider A/bh,s as the prototypical seed 
SMBH mass. 

In our model, we use sin k particles to simulate SMBHs, 
following the prescription of iKrumholz et all |2004 ). When 
the following conditions are met, we create a new SMBH in 
the simulation: 

• The stellar density has to be greater than 0.1 H/cc 
(2.4 x 1O 6 M kpc~ 3 ). This ensures that SMBHs form in 
stellar systems. 

• The stellar 3D velocity dispersion has to be greater than 
100 km/s. With this condition we require that the line-of- 
sight velocity dispersion is am ^ 60 km/s, in agreement 
with the observed Mbh — cr relation. 

• The gas density has to be greater than 1 H/cc. With 
this condition we are sure that seed black holes form in the 
nuclear region of star forming disks. 

• No other sink particle is present within 10 kpc. No new 
seed SMBH will be created within 10 kpc from an old SMBH 
residing at the center of a galaxy. 

Each seed SMBH has a fixed mass Mbh = 10 5 M Q and a 
fixed radius r s i n k = 4Ax ~ 2 kpc, where Aa; is the spa- 
tial resolution in physical units. We assume that the SMBH 
mass is homogeneously distributed inside a sphere of radius 
r s ink, and we add this density distribution to the total mass 
density when solving the Poisson equation. Like dark matter 
and star particles, sink particles are advanced in time by in- 
terpolating the gravitational force back to the sink position 
using the inverse CIC scheme. 

A key ingredient in our simulations is AGN feed- 
back from SMBHs. On the theoretical side, it is thought 
that AGN feedback energy could be transported through 
the gas by seve r al pr ocesses, e.g. radiative feedback 



JCiotti fe Ostriker 2001 ) or strong shocks ( Bcgclman 
12004 iDubois et al.l T20liT iGaspari et al.l 1201 il l. Additional 
feedback may come from production of cos mic rays 
jBriiggen et al.l 120021 ; IChandran fc Raseral 120071 ). On the 
theoretical side, all these processes are very complex to 
model self-consistently and therefore are very challenging to 
be reproduced in hydrodynamical simulations. Additional 
difficulties include the very high spatial resolution required 
to study these processes in a detailed way, although it is not 
only a matter of resolution but also of proper physical mod- 
eling. Even in state-of-the-art cosmological hydrodynamical 
simulations it is not possible to model these processes in a 
self-consistent way, so that a phenomenological treatment of 
the problem must be chosen. 

Mass accretion onto SMBHs and AGN feedback are im- 
pleme nted using a modified version of the lBooth fe Schavel 
l|2009h model which was originally developed for SPH simu- 
lations. This model is a modified version of that proposed by 
ISpringel et all (|2005l V We compute the mass accretion rate 
onto each SMBH using a modified Bondi-Hoyle formula 



Q^boost 



47rG 2 Mg H P 
(cf + u 2 ) 3 / 2 



(3) 



where p, c s and u are the average gas density, sound speed 
and relative velocity with in the sink radius, all computed 
following the approach of IKrumholz et ail (|2004|) . The pa- 
rameter Qboost was introduced by Springel et al.l (|2005l ) to 



account for unresolved multiphase turbulence in the SMBH 
environment, and it's value was originally chosen a s con- 
stant. Numerical studies by iBooth fc Schavel {2009) show 
that Qboost should instead be considered as a function of 
the local density: it should be close to unity in low density 
regions, it should increase in high density regions, in order 
to match the subgrid model used for the unresolved turbu- 
lence in the disks. As proposed bv lBooth fc Schavel {2009), 
we adopt 



Qboost = \~^~J if nH > n * = H/cc, 
Qboost = 1 otherwise. 



(4) 



We stress that the choice for this particular form of Qboost 
is strictly dependent on the chosen EOS for the gas. 

Equation [3] does not provide any upper limit for the 
accretion rate, however it cannot exceed the Eddington limit 



M E 



47rGM B Hra p 

£ r (TTC 



with 



0.1. 



(5) 



where e r is is the efficiency at which accreting gas rest mass 
energy is converted into radiation. To enforce this upper 
limit we always set the accretion rate to 



M acc = min(M B H, Med) 



(6) 



At each time step, a total gas mass of M aC cAt is removed 
from all cells within the sink radius, with the same weight- 
ing scheme as the on e used to define average quantities 
l|Krumholz et al.ll2004l ). In order to prevent the gas density 
from vanishing or becoming negative, we remove no more 
than 50% of the gas at each time step. 

In our simulations, AGN feedback is implemented 
in a way that al l ows s elf-regulated SMBH mass 
growth (ISiiacki et al, I 120071 : ICattaneo fc Tevssierj 120071 : 
IBooth fc Schavell2009l ), by injecting thermal energy directly 
into the gas surrounding the black hole. If the black hole 
mass is too low, the amount of energy injected in the gas 
will not be able to heat it up, so it will remain cold and will 
accrete onto the SMBH at rates close to Eddington limit. 
In this case the SMBH growth can proceed exponentially 
with time, with e-folding time scale equal to the Salpeter 
time 



is = 



e r QTC 
47rGm B 



45Myr 



(7) 



Feedback will become more efficient as long as the SMBH 
continues to grow, until the injected energy is enough to 
unbind the gas surrounding the black hole. In this regime 
accretion will become Bondi-Hoyle limited and mass growth 
will proceed more slowly. Figure [T] shows the time evolution 
of the accretion rate of the most massive SMBH in the sim- 
ulated region. Short bursts of Eddington limited accretion 
are followed by long quiescent epochs of Bondi-Hoyle limited 
accretion, self-regulated by SMBH feedback. 

At each time step we compute the thermal energy in- 
jected in the gas surrounding each black hole as 



AE = e c e r M aC oC At. 



(8) 



where e c is the coupling efficiency, i.e. the fraction of radiated 
energy that is coupled with the surrounding gas. The correct 
value for e c can be set requiring the simulations to reproduce 
the observed Mbh — a relation; we use the fiducial value 
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Figure 1. Time evolution of the SMBH mass in red, and of its 
instantaneous accretion rate in black. 



e c ~ 0.15 (|Booth fc Schavd 120091 ). The energy AE is not 
immediately injected in the gas, but its accumulated and 
stored in a new variable .Eagn, so that we can avoid the gas 
instantly radiating away this energy via atomic line cooling. 
We release this energy within the sink radius when 



-Eagn > 



(9) 



where m gas is the gas mass within the sink radius and T m i n 
is the minimum feedback temperature. T m i n should be cho- 
sen to be at least 10 7 K, the temperature above which line 
cooling is not very efficient, so that the resulting feedback 
is independent of the value of T m i n . In our simulations we 
adopt the fiducial value T m i n = 10 r K. 

This model allows us to phenomenologically reproduce 
the two extreme regimes in which AGN feedback is thought 
to work (|Siiacki et al.ll2007r ): 

(i) " Quasar mode" : in the case of cold dense gas accretion, 
more energy is required to reach the energy threshold. A 
large amount of energy is accumulated and released in a 
burst when the energy threshold is reached. 

(ii) "Radio mode": in the case of hot diffuse gas accre- 
tion, less energy is required to reach the threshold. Energy 
is injected in the gas in a quasi-continuous fashion. 

In the next section we will show how the presence or 
lack of SMBHs and their feedback can heavily influence the 
structure of massive galaxies in clusters. 



3 RESULTS 

In this section, we compare the properties of the BCGs in 
our two simulations to show the strong differences between 
the two models. We also compare our results with the two 
observational sam ples of early-type galaxie s at z ~ and 
2 ~ 1 analysed by Ivan der Wei et al.l |20 0cf). We stress that 
BCGs are particular kinds of early-type galaxies, so we make 



an additional compari son with the BCG sample analysed 
bv lBrough et al.l (|201ll ) to show that our results are robust. 
A final comparison is made with th e cosmological hyd ro- 
dynamical simulations perform ed by | Naab et all J2009I ) (e 



feblCh (the 



massive early- type galaxy) and lFeldmann et al 
central galaxy of a group). 

3.1 Identification of the BCGs 



We identified gal axies in our simula tions with the Adapta- 
HOP algorithm dAubert et all 120041). using the version im- 



plemented and tested bv lTweed et al.l (2009). This version 
of the algorithm allows to use the Most Massive Substruc- 
ture Method (MSM) to identify halos as well as their sub- 
structures, along with their centers, virial radii and masses. 
Galaxy centers were identified by running the same halo 
finder using the star particles. In each of the analysed snap- 
shots, we defined the BCG as the object with the largest 
stellar mass; we verified that such a definition of the BCG 
also selects the galaxy closest to the cluster halo center. 

In both simulations, the BCGs are surrounded by fairly 
smooth stellar halos extending to beyond 50 kpc. Similar 
stellar halos around BCGs have been found in n umerical 
simulations by other authors (|Puchwein et alfeoid ). and are 
similar to the intracluster light observed in galaxy clusters 
jLin et all 120031; iGonzalez et all 120071 ; iGiodini et all |2009| ; 
iRudick et al.ll2010l ). There is no clear separation between 
the BCG and the intracluster stars, so just for accounting 
and comparison purposes, a criterion should be chosen to de- 
cide how to separate the central galaxy from the stellar halo. 
We consider as part of the BCGs all the regions where the 
stellar density is larger than 2.5 x 10 6 Mq /kpc 3 . This choice 
is justified by the fact that lowering the density threshold 
by a factor of two does not influence the values of the mea- 
sured stellar masses by more than 5%. We also checked that 
the mass estimates we obtain using this criterion are robust 
by comparing them with estimates obtained using two addi- 
tional methods: (i) we substituted the 3D density threshold 
with an equivalent surface densi ty threshold ; (ii) w e use a 
criterion similar to that used bv lNaab et alj (|2009l ). where 
the BCG is identified as the set of star particles enclosed in 
a spherical region of fixed physical radius. The variations in 
stellar mass were always less than 5%, therefore we consider 
the masses obtained with the stellar density criterion as our 
fiducial values. 

We find that the BCG mass at z = is M star = 3.484 x 
10 12 M Q in the AGN-OFF run and M star = 2.243 x 10 11 M Q 
in the AGN-ON run. At z = 1 we find M star = 1.637 x 10 12 
M Q in the AGN-OFF run and M star = 1.697 x 10 11 M in 
the AGN-ON run. At both redshifts the total stellar mass of 
the BCG in the AGN-OFF simulations are almost 10 times 
larger. This fact implies that star formation quenching due 
to AGN feedback processes is very efficient in our model and 
has significantly changed the star formation history already 
by a redshift z = 1. 

In Figure [2] we show the comparison of the stellar-vs- 
halo mass relation for the BCGs in our AGN-ON/OFF 
simulations with t he predictions of abundance matching 
|Moster et al.ll2010l ) (black solid line). The BCG mass pre- 
dicted by the AGN-OFF model deviates more than 3a from 
the average relation obtained through abundance matching. 
Without feedback, the BCG forms stars at an efficiency of 
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Figure 2. Comparison of the stellar-vs-halo mass relation in 4 
early-type galaxies from different cosmological simulations (filled 
and empty black dots). The red dotted line represents the rela- 
tion expected for a 20% star formation efficiency from the uni- 
versal baryon fraction. The solid black lin e is the prediction from 
abundance matching (Moster et al. 2010). The grey shaded ar- 
eas represent the la, 2a and 3<r scatter bars around the average 
relation. 

20% from the cosmic available baryon fraction, thus it lies 
close to the dotted red line. On the contrary, the AGN-ON 
model prediction is remarkably close to the average relation, 
and well inside the la scatter bars. For comparison with re- 
cent simulations of massive early-type galaxies that do not 
includ e AG N feedback, we s how t he results of iNaab et al.l 
(|2009h and iFeldmann et all (|2010l ). The mass of the cen- 
tral group galaxy simulated by these authors are both ~ 3a 
above the average relation; also in this case the result is close 
to the simplified model with 20% star formation efficiency. 
This suggests that AGN feedback could well be the mecha- 
nism that reduces the star formation efficiencies in massive 
galaxies. 

3.2 Star formation rates 

The star formation history of the BCG galaxy in our sim- 
ulations can be used to quantify the efficiency of the star 
formation quenching produced by our AGN feedback model. 
The difference between the star formation history between 
the AGN-ON and AGN-OFF simulations is striking. Fig. [3] 
shows the rate at which stars found in the BCG at z = 
have been formed as a function of time. At the beginning of 
the simulation (age of the Universe < 1 Gyr), when SMBHs 
have not formed yet the star formation rates in the AGN- 
ON (blue line) and AGN-OFF (red line) galaxies are al- 
most indistinguishable. Later on AGN feedback starts to 
play its role in quenching star formation: the difference be- 
tween the star formation rate in the two models is almost 



Figure 3. Star formation rate in the BCG as a function of time. 

a factor 10, with t he di fference increasing towards z = 0. 
Irlirschmann et all |201ll) compare star formation rates of a 
wide variety of semi-analytical models and numerical hydro- 
dynamical zoom simulations: they also find that the differ- 
ence between the star formation rates in models including 
AGN and models do not including this effect is striking, 
however smaller t han in our AGN-ON run. A similar result 
has been found bv lvan de Voort et alj (|201ll ). A comparison 
with these results suggests that our AGN feedback model is 
particularly efficient in quenching star formation. It should 
be stressed that the such high feedback efficiency is expected 
because the resolution we reach is sufficient to resolve mas- 
sive galaxies where SMBH reside already at z > 6, thus AGN 
feedback starts playing its role in quenching star formation 
at very early times. It is interesting to note how in both 
the simulations the star formation rate peaks before red- 
shift z = 1 and decreases toward 2 = 0. The fact that star 
formation is extremely quenched between z = 1 and z = 
in the AGN-ON model is also in agreement with the results 
of rec ent semi-analytical models of BCGs (|Lucia Sz Blaizot] 
2007). Some starburst events are evident as peaks in the star 
formation rates at different times, while periods of strong 
AGN activity appear as sudden decreases of the star for- 
mation rate. This simple test shows how our simple AGN 
feedback model is indeed able to quench star formation, de- 
creasing the stellar mass observed in the BCG. 

3.3 Ellipticity and velocity distribution of stars in 
the BCGs 

We constructed mock images of the central region of the 
simulated galaxy clusters to detect further differences be- 
tween the two models. Figure [4] shows a set of such images 
at z = 1 and z — 0, where the color of each pixel represents 
the flux in the g' band. We have chosen the same spatial 
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BCG properties 



Simulation 


Redshift 


M star [M ] 


i? eff [kpc] 


cr cff [km/s] 


e 


v/<T 


AGN-ON 


z 


= 


2.243 


x 10 11 


10.286 


292.96 


[0.123,0.155] 


0.08 


AGN-OFF 


z 


= 


3.484 


x 10 12 


6.858 


652.83 


[0.257,0.498] 


0.65 


AGN-ON 


z 


= 1 


1.697 


x 10 11 


8.254 


277.93 


[0.201,0.424] 


0.16 


AGN-OFF 


z 


= 1 


1.637 


x 10 12 


4.784 


525.32 


[0.120,0.721] 


0.82 



Table 2. Properties of the BCGs. (1) Type of simulation. (2) Redshift z. (3) Stellar mass M B tmt- (4) Half-mass radius R c g. (5) Velocity 
dispersion within R c ff, <T e ff . (6) Ellipticity e at R e ff. (7) Average v/a within R c g. 



and color scales in all the pictures to show the remarkable 
difference in sizes between the AGN-ON and OFF simula- 
tions. The AGN-OFF BCG is much more extended than 
the AGN-ON one, at z = 1 as well as at z — 0, i.e. the 
global size of the galaxies scales with their mass. However, 
the half-light isophotes (black contours) show that the stel- 
lar light is more concentrated in the AGN-OFF BCG than 
in the AGN-ON one. This implies that the surface density 
profile of the AGN-ON galaxy is shallow when compared to 
the AGN-OFF one; this is indeed the case as we will show 
in E31 

To estimate the range of possible ellipticities an ob- 
server would see, we view the simulations from different an- 
gles and make the same surface brightness images. We then 
fit the half-light isophotes with ellipses and estimate their 
ellipticity e = 1 — b/a, where b and a are the semi- minor and 
semi-major axis of the ellipses. In this way we give the range 
of possible ellipticities for our galaxies. At z — the AGN- 
ON BCG has e = [0.123,0.155], while the AGN-OFF BCG 
has e — [0.257,0.498]; both ranges are consistent with val- 
ues typical of spheroidal systems, but the AGN-ON galaxy 
is slightly more spherical. At z = 1 we have e = [0.201, 0.424] 
for the AGN-ON BCG, and e = [0.120, 0.721] for the AGN- 
OFF BCG; at this redshift, the AGN-ON galaxy is slightly 
more spheroidal than at z — 0, while the ellipticity of the 
AGN-OFF galaxy is more typical of flattened ellipsoidal sys- 
tems. 

Figure [5] shows the average line of sight velocity maps 
for the stars in the central region of the cluster at z — 
and z = l. The AGN-OFF BCG hosts a rotating disk with 
a rotational velocity of a few hundreds km/s, both at z = 
1 and z = 0. Interestingly, in the AGN-ON case there is 
no evidence for rotation at z = 0, while at z = 1 there 
is evidence for very slow rotation. These maps should be 
compared with those in Figure [6] that show the line of sight 
velocity dispersion of the stars in the central region. The 
AGN-OFF velocity dispersion at z = is highly peaked 
in the center, where it is larger than ~ 750 km/s, and it 
decreases rapidly towards the external parts of the BCG. On 
the contrary, the AGN-ON velocity dispersion at z = does 
not peak in the center, but it increases towards the outskirts 
of the BCG (we will discuss this topic in greater detail in 
the next subsection). At z = 1 the velocity dispersion is 
centrally peaked in both the simulations, but the values are 
at least a factor of two larger in the AGN-OFF case than in 
the AGN-ON case. If we rotate the galaxies edge-on and we 
account only for the regions within the half-light isophotes 



at z = 0, we measure a mean v/a ~ 0.65 for the AGN-OFF 
BCG, and v/a ~ 0.08 for the AGN-ON BCG. At z = 1 we 
find v/a ~ 0.82 for the AGN-OFF BCG, and v/a ~ 0.16 for 
the AGN-ON BCG. 

It is interesting to compare the (e,v/a) values with 
those observed in early-type galaxies, despite the fact 
that BCGs are a parti cular s ub-category of t h ese o bjects. 
ICappellari et al.l \200H \ and lEmsellem et alJ (|2007l ) used 
2D spectroscopy measurements from the SAURON survey 
l|Bacon et alj|200ll ; IdeZeeuw et al.ll2002h and analysed the 
orbital structure of 48 nearby SO and E galaxies. They divide 
the galaxies in their sample into two groups: slow rotators 
with e < 0.3, v/a < 0.2 and low specific angular momentum, 
and fast rotators with e < 0.7, v/a > 0.2 and a high specific 
angular momentum. Including the effects of AGN feedback 
in our simulations moved our BCG from the class of fast 
rotating early-type galaxies to that of the slow rotators. 

The SAURON galaxies are not chosen to be central 
cluster galaxies, so we can also directly compare our sim- 
ulations with the central galaxy in the Virgo Cluster, M87 
(NGC 4486), that is also a slow rotator; according to 
ICappellari et all (|2007T ). M87 has (e, v/a) = (0.04, 0.02), val- 
ues that are consistent with the AGN-ON BCG simulation. 
This close match with M87 is interesting, given that in our 
simulations we only selected the Virgo-like halo based on its 
virial mass and merger history. 



3.4 Rotational velocity, velocity dispersion and 
stellar mass surface density profiles 

Here we analyse the circular velocity and velocity dispersion 
profiles of the two galaxies at z — (Figure [7J left panel) 
to further explore the structural and kinematic properties 
of our BCG simulations. The rotational velocity (filled sym- 
bols) of the AGN-ON BCG within R eB (blue dashed line) is 
never larger than 25 km/s and its value is slightly growing 
with the distance from the center. The line of sight veloc- 
ity dispersion profile is nearly constant within the BCG, 
250 < a < 300 km/s within R c s. Within this characteristic 
radius the galaxy is always dominated by random motions 
and has a very little rotation. 

The AGN-OFF BCG has very different properties. Its 
velocity dispersion profile is maximum in the central region 
and decreases with radius; the rotational velocity increases 
with distance from the center, peaking just before R e g is 
reached, then decreases again. Within R B ff the AGN-OFF 
galaxy is stil dominated by random motions, however, the 
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Figure 4. Fluxes in the g' band of the BCG projected along the z axis of the periodic box. We filtered out all pixels with flux < 10 17 
erg/s/cm 2 . Top-left: 2 = 0, AGN-OFF. Top-right: 2 = 0, AGN-ON. Bottom-left: 2 = 1, AGN-OFF. Bottom-right: 2 = 1, AGN-ON. The 
black contours represent the half-light isophotes. 



rotational velocity is a significant fraction of the velocity 
dispersion (> 0.45). Whilst rota ting disk-like structur es in 
fast rotators have been observed |Krainovic et al.ll200&t ). the 
peak u rot value of our AGN-OFF simulation is a factor ~ 
2 — 3 larger than the typical values measured for the fastest 
rotators. 

The veloc ity dispersion profiles measured for M87 

jMurphv et al l l201ll ) are comparable with those of our 
AGN-ON BCG, at least beyond the spatial resolution limit. 
In SAURON the kinematic properties of M87 h ave been 
measured down to ~ 0.1 kpc |Emsellem et al.ll2004h . but our 



spatial resolution limit is much larger: our AGN-ON model 
is not able to reproduce the velocity dispersion profile ob- 
served in M87 for R < 1 kpc. However, when we consider the 
resolved region, R > 1 kpc, the veloci ty dispersion profile of 
our AGN-ON mode l and that of M87 (|Emsellem et al.ll2004l : 
iMurphv et aHl201lD are rather similar. 

On the right panel of Figure [7| we show the stellar mass 
surface density profiles at z = for our two BCGs. We find 
that these profiles can be fitted by Sersic profiles for R > 10 
kpc for the AGN-ON BCG (blue line, Sersic index n = 10) 
and R > 4 kpc for the AGN-OFF BCG (red line, Sersic index 
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Figure 5. Average velocity along the line of sight. The projection is along the 2 axis of the periodic box. Top-left: 2 = 0, AGN-OFF. 
Top-right: 2 = 0, AGN-ON. Bottom-left: 2 = 1, AGN-OFF. Bottom-right: 2 = 1, AGN-ON. 



n — 10). We find that there is almost an order of magnitude 
difference between the stellar mass surface density of our 
two models. It is also interesting to compare the measured 
profiles with the Sersic functions at 1 kpc < R < 10 kpc: in 
the AGN-OFF case we find that there is a mass excess with 
respect to the Sersic function and in the AGN-ON case we 
find a mass deficiency. On the observation al side, similar re- 
sults were describe d by several autho rs (lKormendvl | 19991 : 
Guillen et all l2000l; lLaine et alj 120031: iTruiillo et aljpool 
Lauer et all 120051: ICdte et al.l [20071; iKormendv et all 120091 



Graham 2011 V Kormendv et al.f i 20091 ) analyse the surface 
brightness profiles of elliptical and spheroidal galaxies in the 



Virgo cluster and show that it is possible to fit them with 
single light profiles to large radii. The fits are generally not 
satisfactory in the central regions; in particular, they find 
that most luminous elliptical galaxies have either light ex- 
cesses or deficiencies at there centres. Light deficiencies are 
preferably found in the most luminous objects like NGC 
4472 and M87 and the regions where they are found extend 
from the center to a distance < 10 kpc. If we assume that 
our surface density profiles can be mapped directly to sur- 
face brightness profiles, the observations of luminous galax- 
ies are a much closer match to the AGN-ON simulation. In 
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Figure 6. Velocity dispersion along the line of sight. The projection is along the z axis of the periodic box. Top-left: 2 = 0, AGN-OFF. 
Top-right: 2 = 0, AGN-ON. Bottom-left: 2 = 1, AGN-OFF. Bottom-right: 2 = 1, AGN-ON. 



the next subsection we discuss the core formation process in 
more detail. 



3.5 Formation of the stellar core 

The stellar surface density profile of the AGN-ON BCG 
within the inner 10 kpc has constant density core. In the 
left panel of Figure [8] we show the time evolution of the 
stellar surface density profile of the AGN-ON BCG, from 
2 = 2 to z — 0. From this plot we can see that early on 
the stars are centrally concentrated, but a significant core 
starts developing at a redshift z=l. The core seems to be 



in place by a redshift z=0.6, or about 3 Gyrs, which is a 
cluster crossing timescale, but about 100 crossing times at 
the half light radius of 10 kpc. The core in the AGN-ON 
BCG doesn't form naturally, the AGN-OFF simulation has 
no such effect. Something has caused stars to move away 
from the inner region and this requires a significant amount 
of energy. 

iNaab et al.l (|2009l ) showed that repeated minor dry 
mergers lead to large increases in the size of massive elliptical 
galaxies and to decreases of the central stellar density con- 
centration. The cores of accreted galaxies act as perturbers 
and modify the stellar mass distribution through dynami- 
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Figure 7. Left panel: stellar velocity dispersion and rotational velocity of the BCG in our two models. The black solid lines are for the 
AGN-ON model and the black dashed lines are for the AGN-OFF model. The dotted lines indicate the effective radii of the AGN-ON 
(blue) and AGN-OFF (red) BCGs. Right panel: stellar mass surface density profile of the BCG in our two models, compared with Sersic 
profiles (see text for details). The black solid lines are for tha AGN-ON model and the black dashed lines are for the AGN-OFF model. 
In both panels the grey shaded area shows the unresolved region of our simulations. 




Figure 8. Secular evolution of the stellar surface density profile towards a cored profile (left panel). The right panel shows the cumulative 
SMBH mass, showing that the late time evolution is governed by a succession of dry SMBH mergers towards the central galaxy, while, 
at early time, the total mass in SMBH is still growing by gas accretion. 



cal friction. We do not see this effect in our AGN-OFF run 
because completely dry mergers are rare due to the large 
reservoirs of gas residing in galaxies. In the AGN-ON case a 
significant fraction of the ga s is removed from galaxies and 
the mechanism proposed bv lNaab et al.l (|2009l ) can be effi- 



cient. As we will show in Subsection 13. 6 1 we observe that the 
AGN-ON BCG is much mor e extended than th e AGN-OFF 
BCG, in agreement with the lNaab et al.l (|2009h picture, but 
it is challenging to produce an extended flat core like the 
one observed in our AGN-ON model only through repeated 
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dry minor mergers. SMBHs and AGN feedback provide ad- 
ditional processes able to contribute to the core formation 
mechanism. 

SMBH binaries are expected to form naturally in 
the hierarchy of mergers that lead to massive ellipti- 
cals and they provide several mechanisms that can pro- 
duce cores. The most favoure d model is SMBH scouring 
l|Milosavlievic fc Merrittll2003h . At the galaxy centres the 
black holes form binary pairs that decay as they transfer 
energy to stars via three body encounters, thus they are 
able to eject stellar material from the central regions and 
form a core. Numerical experiments suggest that the cumu- 
lative effect of multiple SMBH dry mergers is able to remove 
a stellar mass that is ~ 2 — 4 times the final SMBH mass 
l|Merritt et al.l l2007). This process is important from about 
100 to 1 parsecs, so it is unresolved in our simulations which 
have 1 kpc force softening. At the softening length, the en- 
closed stellar mass is larger than the BH masses therefore 
a binary BH system can't form. However there is an ad- 
dition process that we do resolve: During the mergers the 
black holes would sink to the very central region due to 
dynamical fri c tion. The numerical experiments presented in 
IGoerdt et all l|201Ch show that the energy transferred from 
the sinking SMBHs to stars contributes to the formation of 
cored profiles; the typical mass deficits are found to have a 
similar magnitude as the SMBH mass. 

The right panel of Figure[5]shows the cumulative SMBH 
mass from the center out to virial radius of the cluster at 
different output times. Each discontinuity in these profiles is 
associated to the presence of a SMBH. At early times the to- 
tal mass in SMBHs is increasing mainly due to gas accretion, 
but at later times the growth is dominated by a succession 
of SMBH mergers. The mass of the central S MBH at z = 
is Mb h = 4.2 x 10 9 Mo The simulations of IGoerdt et all 
|2010l ) and iMerritt" et al. (2007) used idealised equilibrium 
models to study the effects of sinking massive binary objec ts. 
According to the model developed in lGoerdt et al.( l|201Cf ). a 
sinking SMBH of mass Mbh = 4.2 x 10 9 M Q would produce 
a core of size R CO ic ~ 3 kpc. 

It should be stressed that these dynamical effects of 
sinking SMBHs act on all components of the mass distri- 
bution, although how stars, gas and dark matter respond 
can be quite different. Given this caveat, we note that the 
predictions for the mass deficit in such models is about one 
third lower than the mass deficit within the inner 10 kpc 
of our AGN-ON BCG, M de f = 3.04 x 10 10 M , measured 
as the difference between the mass predicted by the Sersic 
profile in Fig. [7] and the actual mass enclosed in the same 
region. However, recent N-body experiments performed by 
Kulk arni fc Loeq (|201 lh showed that when multiple (three 
or more) SMBHs are present, the core formation process is 
much more efficient: mass deficits in this case can be more 
than 5 times the total SMBH mass. These results bring the 
mass deficit obtained in pure N-body simulations with per- 
turbing SMBHs much closer to the mass deficit observed in 
our AGN-ON run. 

Additional energy for the formation of the core can be 
provided by strong AGN driven outflows at z < 1 that mod- 
ify the local gravitational potential and may cause expan- 
sion of both the dark and stellar mass distribution. Sim- 
ilar processes have been observed in numerical cosmolog- 
ical simulations of dwarf galaxies in which gas outflows 



are generated by supe rnovae feedback (|Navarro et al.lll996l ; 



iGovernato et al 



2010l). In the N-body simulations presented 



in lNavarro et al ( 19961 ). the mass outflows are simulated by 
growing and rapidly removing an idealised potential from 
the centre of an equilibrium realisation of a dark matter halo; 
the effect is the formation of a core in the density profile. 
In this model the efficiency of the process is tx M^^R^^ 2 , 
where Mdi sc is the mass of the disc and -Rdisc is its scale 
radius. This model is substantially different from ours, since 
we have multiple epochs of AGN driven gas outflows. How- 
ever, these gravitational potential fluctuations should act 
in a similar way and could contribute to creating the large 
stellar cores we find. 

The maps in Fig. [9] show the structure of such out- 
flows in our AGN-ON simulation. We show the gas prop- 
erties just after a strong AGN burst at z ~ 0.1: the map 
in the left panel shows the radial velocity of gas with re- 
spect to the cluster centre (identified as the centre of the 
BCG), while the right panel shows a map of the gas entropy 
K = k^T/n 2 ^ 3 . These plots show that in the outer regions 
of the cluster high entropy gas is fastly moving towards the 
center. This corresponds to hot intracluster medium cooling 
and flowing towards the cluster center. On the contrary, the 
inner regions (r <400 kpc) are characterized by convective 
motions pushing gas outwards at velocities ~ 200 km/s, e.g. 
note the ring-like structure observed in the radial velocity 
map. Looking at the entropy map we see that gas within 
400 kpc from the centre is typically low entropy material. 
Observing the regions very close to the cluster centre and 
the BCG, we can speculate on the origin of the convective 
motions observed at higher distances. Close to the BCG a 
significative fraction of the gas is falling towards the centre, 
but it is possible to observe an extended bubble of gas mov- 
ing away from the centre at speed > 400 km/s. This bubble 
has particularly high entropy and corresponds to gas that 
has been heated by a strong AGN burst. This gas will be 
eventually pushed away from the central region and will mix 
with the hot intracluster medium. The net effect of AGN ac- 
tivity is that gas is continuously expelled from the cluster 
inner regions. 

The plot in Figure [TU] highlights the effect of gas out- 
flows in our simulations: we plot the stellar (solid line) and 
gas (dashed line) mass enclosed within 10 physical kpc as a 
function of redshift for both the AGN-ON (blue) and AGN- 
OFF (red) simulations. There is more than an order of mag- 
nitude difference between the central stellar/gas mass in the 
two simulations. The plot shows that in the AGN-OFF case 
gas continuously accretes onto the central region, while the 
stellar mass is constantly increasing due to star formation 
events; the sudden increase of the gas mass are associated 
with galaxy mergers. A different scenario emerges in the 
AGN-ON simulation: the central gas mass is roughly con- 
stant at z > 1, due to the fact that AGN feedback is efficient 
in heating and pushing gas away. In this regime star forma- 
tion is still proceeding but it is strongly quenched. At z < 1 
AGN feedback in the central region becomes extremely ef- 
ficient and ~ 10 11 Mq of gas is expelled from the center 
before z = 0. At < z < 1 the stellar mass inside the core is 
almost constant because there is very little star formation. 
The large gas outflows generate fluctuations in the gravi- 
tational potential that cause the ejection of stellar mass, 
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Figure 9. Left panel: gas radial velocity map after an AGN burst in the AGN-ON simulation. Only the inner part of the cluster is shown. 
Positive values (green to red) of the radial velocity correspond to outflowing material, negative values (blue to violet) are associated to 
inflowing material. Right panel: gas entropy map for the same region. 
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Figure 10. Time evolution of the total stellar and gas mass 
enclosed within 10 physical kpc. 



kpc, which is close to the 10 kpc core observed in our simu- 
lations. 

Finally, we note that the gas mass within the cluster 
core of the AGN-ON simulation slowly decreases over time. 
This slow decrease in gas mass will lead to an adiabatic ex- 
pansion of the total mass distribution, which will also con- 
tribute to the formation of a central core. Furthermore, the 
cooling time of hot gas within the inner 30kpc of the cluster 
centre is about one Gyr, thus we envisage that in its quiet 
mode the AGN can slowly eject the gas that rains down onto 
the centre from the inner cooling flow. 

To further test the efficiency of these processes, ex- 
tensive numerical tests should be performed using ad hoc 
high resolution N-body simulations. We will go into further 
details about how these processes manifest themselves in 
galaxy clusters in a future paper (in preparation). 

As mentioned in the previous section the size of the 
core in our AGN-ON BCG is slightly larger than what 
is expected for typical luminous elliptical galaxies. Also, 
the central surface brig htness profiles of ellipticals in Virgo 
l|Kormendv et al1l2009l ) are seldom as flat as the stellar pro- 
file we observe in our AGN-ON BCG. This may suggest that 
the model we adopt to describe AGN feedback is slightly too 
efficient in injecting energy in the gas and in removing mass 
from regions surrounding the SMBHs. 



possibly preventing the formation of a cusp in the surface 
density profile. 

The efficiency of core formation through such mass 
out flows has not been studied before. Naively applying 
the lNavarro et all i 1 !)!)(>) model under the assumption that 
Miisc equals the ejected gas mass and that i?disc = 10 kpc, it 
is possible to predict the formation of a core of size R corc ~ 6 



3.6 Masses, sizes and velocity dispersions 

We now compare the glob al properties of the galax ies with 
the observational data of Ivan der Wei et ail (|200Sl ). Their 
first sample is composed of > 17000 early- type galaxies at 
redshifts 0.04 < z < 0.08 taken from the SDSS database and 
whose stellar masses, effective radii (recovered from SDSS 
g-band imaging) and average velocity dispersion within the 
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effective radii are known. The second sample is composed 
of 57 cluster galaxies at redshifts 0.8 < z < 1.2, also with 
known stellar masses, effective radii (recovered from zsso- 
band imaging) and average velocity dispersions. Some of the 
most massive early-type galaxies in these samples are BCGs. 
To make our comparison with observations more focused on 
BCGs, we also compare our simulat ed central gala x ies to the 
four BCGs at z < 0.09 analysed bv lBrough et al.l (|201lh . 

For the two simulated BCGs we measured stellar 
masses, M sta r, half-mass radii, .Reflu, and velocity disper- 
sions within the half- mass radius, u (t at z = and 2=1. 
The results are summarized in Table [3711 The difference be- 
tween the properties of our two BCGs is remarkable at z = 
as well as at z — 1. We find that: 

• At z = 1 the AGN-OFF stellar mass is an order of 
magnitude larger than in the AGN-ON case. The difference 
is even larger at 2 = 0. 

• The core of the stellar distribution, quantified by R c g, 
is almost twice as extended in the AGN-ON case than in 
the AGN-OFF case, both at z = 1 and 2 = 0. 

• The average velocity dispersion within R e fi is a factor 
~ 2 larger in the AGN-OFF case at z = as well as at z = 1. 



The results of iNaab et all (|2009T i imply that the prop- 
erties of very massive early-type are determined by repeated 
dry minor mergers that lead to the increase of galaxy sizes 
and masses between z — 1 and z = 0. In our AGN-ON model 
a large fraction of the gas is expelled from galaxies due to 
AGN feedback, increasing the probability of dry mergers. 
This effect is completely lacking in the AGN-OFF run be- 
cause much more gas is present in all galaxies. This may 
partially explain why the AGN-ON BCG is more extended 
than its AGN-OFF counterpart. 

Figure [11] shows the distribution of galaxies at low red- 
shift in the (cr e ff, Red) plane (left) and in the (M s tar, -Reff) 
plane (right). We c ompare the observ ation s with our simu- 
lation s and those of INaab et al.l |2009l 'l and iFeldmann et al.l 
|20ld ) (which do not include SMBHs nor AGN feedback). 
The simulation performed bv lNaab et al.l (|2009T ) reproduces 
quite well the properties of a typical early-type galaxy em- 
bedded in a Milky- Way sized halo at z = 0. The c entra l 
galaxy in the group simulated by IFeldmann et al.l (|201Cj ) 
has a stellar mass and effective radius compatible with sev- 
eral observed early- type galaxies, but a velocity dispersion 
slightly larger than expected. 

The AGN-OFF BCG of this paper is completely dif- 
ferent from any observed low redshift massive early-type 
galaxy in the SDSS sample. Its mass and velocity disper- 
sion are larger than those of any observed early-type galaxy 
by a factor ~ 3 and ~ 2, respectively. On the contrary, the 
properties of the AGN-ON BCG are quite consistent with 
those of observed galaxies. We point out that the proper- 
ties of o ur AGN-ON BCG a re also very similar to the four 
BCGs of lBrough etHI |201ll ). These facts give strong obser- 



1 We assume that the radius that encloses half the mass of stars 
is close to the half light radii. 



vational support to galaxy formation models including AGN 
feedback. 

To further test the effects of including the AGN we 
study the properties of the simulated galaxies at z = 1. 
Figure [12] shows the distribution of the sample of cluster 
galaxies at z ~ 1 in the (a e g,R e fi) plane (left) and in the 
(Mstar, R e fi) plane (right) . For c omp arison, we also show als o 
the results of lNaab et ' all l|2009h and lFeldmann et all < |2010h . 
even if they do not refer to galaxies in clusters. Again, we 
find that the AGN-OFF BCG does not match the observa- 
tions. The properties of the AGN-ON BCG at z = 1 are 
quite close to those of some of the most massive cluster 
galaxies in the sample, but the agreement is worse than at 
low redshift. 



4 SUMMARY AND CONCLUSIONS 

We simulated the formation of a Virgo-sized galaxy cluster 
to study the properties of the massive bright galaxy that 
assembled at its centre. We considered two different mod- 
els: the first one includes standard galaxy formation recipes 
(gas cooling, star formation and supernovae feedback) but 
no AGN feedback, whereas the second one also includes the 
effect of AGN feedback. Our results suggest that the ef- 
fects of AGN feedback in clusters of galaxies is very impor- 
tant for the forma tion a nd evolution of t h e BCG s. Whereas 
INaab et al l (|2009T l and IFeldmann et~afl |20ld ) show that 
AGN feedback is not strictly required to reproduce the 
properties of massive early-type galaxies in the field and in 
groups, we provide evidence that AGN feedback is needed to 
form realistic BCGs in clusters, i.e. in more massive halos. 

We compared the BCGs of our two simulations and 
we found substantial differences: when AGN feedback is 
neglected we obtain an elliptical galaxy whose properties 
(mass, size, kinematic structure, stellar density profile) are 
incompatible with observed elliptical galaxies at z ~ and 
2 ~ 1. It is an extremely massive and fastly rotating galaxy, 
with a stellar cusp in the center. When AGN feedback is 
included the BCG appears to be completely different: it is 
10 times less massive because star formation quenching is 
very efficient; it slowly rotates and its stellar surface density 
profile is cored in the inner 10 kpc. 

AGN feedback results in a stellar-to-halo mass ratio 
is consistent with the predicition of abundance matching 
jMoster et alj|20ldl . A comparison with the massive early- 
type g alaxies in the samples analysed by Ivan der Wei et al.l 
(2008) shows that the mass, the velocity dispersion and the 
effective radius are consistent with those of the most mas- 
sive early-type galaxies observed in the SDSS at z ~ 0, and 
cluster galaxies at 2 ~ 1. We note that a slight decrease 
in the efficiency of AGN feedback would produce a slightly 
larger mass and a lower effective radius at z — 1, bringing 
our simulated galaxy into an even closer agreement with the 
observations. 

The existence of the core in the stellar surface density 
distribution is in agreement with what is observed for the 
most luminous and massive galaxies in the Virgo cluster 
that s how significant mass deficiencies in their central re- 
gions 



Truiillo et al 1 l2004l : iCote et al 1l2007l : lKormendv et al.l 



2009; Graham I2011T ) . We have discussed several mechanisms 
that could contribute to the shaping of the final properties 
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Figure 11. Velocity-size (left panel) and mass-size (right panel) relation of early- type galaxies at rcdshift z = from Sloan data, 
compared to four early- type galaxies from different cosmological simulations. The black dotted lines are contours of the number of 
early-type galaxies per bin in the 0.04 < z < 0.08 sample; going from outside-in we show contours for 5, 10, 30, 100, 200, 300, 400, 500, 
600 galaxies per bin. Each bin has a size Alog(Mstar) = Alog(i? c ff) = Alog(cr c ff) = 0.1. The four BCGs analysed by Brough et al. 2011 
are also shown as black filled squares. 
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Figure 12. Velocity-size (left panel) and mass-size (right panel) relation of clusters elliptical galaxies at redshift 2 = 1 from van der 
Wei (2008), compared to 4 early- type galaxies from different cosmological simulations. 



of the BCG and, especially, to the formation of its core: 
(I) a series of dry mergers that lead to SMBHs sinking to 
the halo center via dynamical friction. This process can 
eject a large fraction of sta rs and dark matter from the 
central regions of the BCG (|Milosavlievic fc Merritdl2003l ; 



iGoerdt et aT1l2010l ). (II) AGN feedback driven gas outflows 
can modify the gravitational potential in the regions close 
to SMBHs; these outflows are impulsive and the 'revirialisa- 
tion' of the inner mat erial can lead to the formation of a core 
l|Navarro et al.|[l996l) . (Ill) The central hot gas slowly cools 
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radiatively, falling onto the SMBH in a convective flow and 
is subsequently ejected impulsively. The slow loss of mass 
from the central region will result in the inner mass distribu- 
tion expanding. The efficiency of each of these mechanisms 
will be explored using idealised numerical experiments in a 
subsequent study. 

Observations show that low mass early-type galaxies 
typically have cusps in their surface brightness profiles, 
while high mass ear ly-type galaxies preferably have cen- 
trally cored profiles touiillo et al l 12004 ICote et all 120071 ; 
iKormendv et al] 120091 ) . We find that neglecting the pres- 
ence of SMBHs and AGN feedback produces a cusp, while 
including these effects produces a core. These considerations 
suggests that there may be a close connection between the 
mass dichotomy in early-type galaxies and the presence of 
SMBHs. In high mass early-type galaxies the efficiency of 
the processes that lead to a core formation are expected to 
be higher than in lower mass early-type galaxies, thus lower 
mass galaxies may retain the cusps in the distribution of 
their stars. 
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